ccccccccc
c  CSHOT2 builds traces given the ray data calculated by CSHOT1.
c  There are 2 modes of operation:
c  (1) Generate shot gathers
c  (2) Generate receiver gathers

c  For case (1), you specify the record numbers of the shots to build.
c  Also, the traces to generate within each shot are given.  You can
c  create a common offset data set by choosing to build only 1 trace/shot
c  and by being careful with the geometry cards in CSHOT1 (ie, you have to make
c  sure that the trace in question is the same distance from the shot in all
c  records).

c  In case (2) you specify a range of shot records to search for the common
c  receiver data.  You also specify the range of receivers for which gathers
c  are to be built.
ccccccccc


cc Variable List:

c  AMP      Wavelet is caled by this amplitude factor (calculated   
c           by CSHOT1)
c  CAUSTC   TRUE if the ray has passed through a caustic (wavelet  
c           will be phase shifted)
c  DATIN    Unit number of input data (generated by CSHOT1)
c  DATOUT   Unit number of output data (traces)
c  DELTAT   Sample rate used in calculating wavelet
c  EOF      TRUE if the end of file of the input data has been reached
c  FLL,FHL,FLR,FHR  Amplitude spectrum of wavelet
c  FNYQ     Nyquist frequency
c  FUDGE    Constant scale factor applied to all amplitudes
c  HEAD     TRUE if this is a head wave event
c  I        Counter
c  IBACK    Number of records to backspace in input data file to get
c           back to start of current shot (used when paneling the output)
c  IDUMMY   Dummy integer variable
c  INFIL    Name of the input data file (created by CSHOT1)
c  IRECD    Record number
c  IRECD1   First (shot) record in sort
c  IRECDN   Last (shot) record in sort
c  IRECV    Station number of the receiver
c  IRPREV   Previous record number
c  IROUT    Loop variable over receiver sorts
c  IROUT1   First receiver sort to output
c  IROUTN   Last receiver sort to output
c  ITR1     First trace to build
c  ITRN     Last trace to build
c  ISOUT    Current output record number
c  ITRACE   Trace number
c  ITRMAX   Maximum trace to output in this panel of the shot record
c  ITRMIN   Minimum trace to output in this panel of the shot record
c  J,K,L    Counters
c  JOB      Job descriptor (s, r, or o)
c  KEY      Either the recievr number or the trace number, depending
c           on the sort option
c  MAXTPT   MAximum number of samples/trace
c  MAXTR    Maximum number of traces/panel - for shot records
c  MAXWAV   MAximum number of samples in the wavelet (at input sample rate)
c  NMOD2    Number samples in wavelet modulo 2
c  NPTBY2   NPTWAV / 2
c  NPTWAV   Number of time samples in wavelet (at input sample rate)
c  NREC     Number of traces/shot in input data set
c  NRSWAV   Number of points in the resampled wavelet
c  NSHOTS   Number of shots in the input data set
c  NTPTS    Number of samples/trace
c  NTRACE   Number of traces per output panel
c  NWAVMX   Maximum number of points in resampled wavelet
c  OFFSRT   TRUE for common offset (ie, common trace) sort
c  OUTFIL   Name of the output file containing the traces
c  PARIN    Unit number of input parameter file PARAM2
c  PHASE    Phase due to post critical reflections
c  PI       3.1415927E+000
c  RECSRT   TRUE for a receiver sort
c  RESAMP   Factor by which to resample wavelet (should be even)
c  SHTSRT   TRUE for shot sort
c  SRATE    Sample rate in seconds
c  TARGET   Receiver number or trace number for sort
c  TIME     Arrival time of wavelet
c  TMAX     Trace end time
c  TRACE(,) Array containing traces 
c  WAVLEN   Length in seconds of the wavelet
c  WHEAD()  Head wave wavelet
c  WIMAG()  Imaginary part of (zero phase) wavelet
c  WREAL()  Real part of (zero phase) wavelet    10240    5104
c  XDUMMY   Dummy float variable
c  ZDUMMY   Dummy float variable


      INTEGER     DATIN,    DATOUT,    parin,    STDERR

      PARAMETER ( STDERR = 0,
     :            DATIN  = 12,
     :            DATOUT = 13,
     :            parin  = 14)

      INTEGER     MAXTR,        MAXTPT,         MAXWAV

      PARAMETER ( MAXTR  = 500,
     :            MAXTPT =1001,
c    :            MAXWAV = 513)
     :            MAXWAV = 257)

      integer     resamp,  nwavmx
c     resamp should be an EVEN integer OR 1 for no resampling
      parameter ( resamp = 64)
      parameter ( nwavmx = maxwav*resamp )

      REAL       WREAL(2*nwavmx),   WIMAG(nwavmx),
     :           whead(nwavmx),     TRACE(MAXTPT,maxtr)

      REAL       AMP,      deltat,  FUDGE,
     :           FLL,      FHL,     FLR,     FHR,    FNYQ,
     :           PHASE,    PI,      SRATE,   TIME,   TMAX,
     :           WAVLEN,   xdummy,  zdummy

      INTEGER    I,        J,       K,       L,      iback,
     :           isout,    irecd1,  irecdn,
     :           irecd,    irprev,  irecv,   key,
     :           itrace,   itrmin,  itrmax,  ntrace,
     :           nmod2,    NPTWAV,  nptby2,  NREC,   nrswav,
     :           nshots,   NTPTS,   itr1,    itrn,
     :           irout,    irout1,  iroutn,  idummy


      CHARACTER   INFIL*20,   OUTFIL*20,   job*1

      logical     shtsrt,  recsrt,   caustc,   head,
     :            eof,     offsrt

      PARAMETER ( PI = 3.141592653589,
     :            FUDGE = 1.)



      OPEN(UNIT=parin,FILE='param2',STATUS='old',ERR=5)
      GO TO 10
5     WRITE(STDERR,'(A)') 'Can''t open file : ''param2'''
      STOP

10    REWIND parin
c     read job description

      READ(parin,'(A)',err=15) job
      go to 16
15    write(stderr,'(1x,a)')
     :'MAIN: Can''t read job option, file ''param2'''
      stop
16    continue

c     read first and last records to sort on 
      read(parin,*,err=18) irecd1, irecdn
      go to 19
18    write(stderr,'(1x,a,a)')
     :'MAIN: Error reading first, last records for sort,',
     :' file ''param2'''
      stop
19    continue

      if(job.eq.'s'.or.job.eq.'S') THEN
         shtsrt = .true.
c        read first and last traces to build
         read(parin,*,err=20) itr1,  itrn
         go to 21
20       write(stderr,'(1x,a,a)')
     :   'MAIN: Error reading first, last traces to build,',
     :   ' file ''param2'''
         stop
21       continue
      else if(job.eq.'r'.or.job.eq.'R') THEN
         recsrt = .true.
c        read receiver numbers for sort
         read(parin,*,err=23) irout1,  iroutn
         go to 24
23       write(stderr,'(1x,a)')
     :   'MAIN: Error reading receiver sort numbers, file ''param2'''
         stop
24       continue
      else
         write(stderr,'(/1x,a,/33x,a/)')
     :   'MAIN: unrecognized option - use s for shot records,',
     :   'r for receiver sort.'
         stop
      end if


c     get corner frequencies of wavelet
      READ(parin,*,err=30) FLL,FHL,FLR,FHR
      go to 31
30    write(stderr,'(1x,a)')
     :'MAIN: Error reading filter frequencies, file ''param2'''
      stop
31    continue
      IF(FLL.LT.0.) THEN
         WRITE(STDERR,'(1x,A)')'MAIN: negative frequency in spectrum'
         WRITE(STDERR,'(1x,A)')'Specify positive half of spectrum only'
         STOP
      END IF

c     fll<=fhl<=flr<=fhr
      if(fhl.lt.fll) fhl = fll
      if(fhr.lt.flr) fhr = flr
      if(flr.lt.fhl) then
         write(stderr,'(1x,a)')
     :   'MAIN: Error in frequency specification of wavelet.'
         stop
      end if

c     length of wavelet
      READ(parin,*,err=35) WAVLEN
      go to 36
35    write(stderr,'(1x,a)')
     :'MAIN: Error reading wavelet length, file ''param2'''
      stop
36    continue

c     sample rate in seconds
      READ(parin,*,err=40) SRATE
      go to 41
40    write(stderr,'(1x,a)')
     :'MAIN: Error reading sample rate, file ''param2'''
      stop
41    continue
c     resampling rate for wavelet
      deltat = srate / resamp

c     Check frequencies less than Nyquist.
      FNYQ = 1. / ( 2. * SRATE)
      IF(FHR.GT.FNYQ) THEN
         WRITE(STDERR,'(1x,A)') 'MAIN : frequency beyond Nyquist.'
         STOP
      END IF

c     number of points in wavelet
      nptwav = wavlen / srate + 1
      nmod2  = mod(nptwav,2)
      if(nmod2.eq.0) then
c        make number of points in wavelet an odd number
         nptwav = nptwav + 1
      end if
      if(resamp.eq.1) then
         nrswav = nptwav
      else
         nrswav = nptwav * resamp + 1
      end if
      nptby2 = ( nptwav - 1 ) / 2

      IF(nrswav.GT.nwavmx) THEN
         WRITE(STDERR,'(1x,A)') 'MAIN : too many points in wavelet.'
         STOP
      END IF


c     trace length
      READ(parin,*,err=50) TMAX
      go to 55
50    write(stderr,'(1x,a)')
     :'MAIN: Error reading trace length.'
      stop
55    continue

      NTPTS = 1 + ( TMAX + SRATE / 2. ) / SRATE
      IF(NTPTS.GT.MAXTPT) THEN
         WRITE(STDERR,'(A)') 'MAIN : too many points / trace.'
         STOP
      END IF

      IF(NPTWAV.GT.NTPTS) THEN
         WRITE(STDERR,'(A)') 'MAIN : wavelet length > trace length.'
         STOP
      END IF

c     Read the names of the input and output data files
      READ(parin,'(A)',err=70) INFIL
      go to 75
70    write(stderr,'(1x,a)')
     :'MAIN: Error reading input file name in file ''param2'''
      stop
75    continue

      READ(parin,'(A)',err=80) OUTFIL
      go to 85
80    write(stderr,'(1x,a)')
     :'MAIN: Error reading output file name in file ''param2'''
      stop
85    continue

      OPEN(UNIT=DATIN,FILE=infil,STATUS='old',
     :form='UNFORMATTED',ERR=140)
c    :ERR=140)
      GO TO 145
140   WRITE(STDERR,'(1x,A)') 'Error opening input data file.'
      STOP
145   REWIND DATIN

c     OPEN(UNIT=DATOUT,FILE=OUTFIL,ERR=150)
      OPEN(UNIT=DATOUT,FORM='UNFORMATTED',FILE=OUTFIL,
     :ERR=150)
      GO TO 155
150   WRITE(STDERR,'(1x,A)') 'Error creating output trace file.'
      STOP
155   REWIND DATOUT




c     Constructing the wavelet.
      CALL WAVLET(FLL,FHL,FLR,FHR,NrsWAV,deltat,
     :nwavmx,PI,WREAL,WIMAG,whead,stderr)


      if(shtsrt) then
c        SHOT GATHERS

c        initialize shot counter
         isout = irecd1
         iback = 0

         nrec = itrn - itr1 + 1
c        number of traces/record = number of receivers/shot
900      itrmin = itr1
         if(nrec.gt.maxtr) then
c           will have to panel the output
            itrmax = itrmin + maxtr - 1
            ntrace = maxtr
         else
            itrmax = itrn
            ntrace = nrec
         end if

c        zero the trace array for this shot
         DO 950 i = 1,  ntrace
            do 925 j = 1, ntpts
               trace(j,i) = 0.
925            continue
950         continue

         eof = .false.
1000     read(datin,end=1050) irecd,itrace,irecv,time,
     :   xdummy,idummy,zdummy,
     :   amp,phase,caustc,head

         if(irecd.lt.isout) go to 1000
         if(irecd.gt.isout) go to 1100
         iback = iback + 1
         if(itrace.lt.itrmin.or.itrace.gt.itrmax) go to 1000

c        trace number within this panel
         itrace = itrace - itrmin + 1

         amp = amp * fudge
c        put the wavelet on the trace
         call addwav(phase,caustc,nptwav,nptby2,time,ntpts,trace,
     :               amp,wimag,wreal,whead,pi,itrace,srate,
     :               deltat,head,resamp,maxtr,maxtpt)

         go to 1000

1050     eof = .true.
         if(irecd.lt.irecd1) then
c           at eof and did not find any records within range
c           close file, stop
            write(stderr,'(1x,a)') 'No shots within specified range.'
            write(stderr,'(1x,a,i4)') 'Max record is ',irecd
            close(unit=datout)
            stop
         end if

1100     continue
c        output the shot panel
         write(stderr,'(1x,a,i4,1x,a,i4,a,i4)')
     :   'Outputting shot ',isout,'traces',itrmin,' - ',itrmax
         DO 1300 I = 1,  ntrace
            write(datout)(trace(j,i),j=1,ntpts)
1300        CONTINUE

         if(itrmax.lt.itrn) then
c           more traces for this shot
            itrmin = itrmin + maxtr
            itrmax = itrmax + maxtr
            if(itrmax.gt.itrn) itrmax = itrn
            ntrace = itrmax - itrmin + 1
c           backspace to beginning of data for this shot
            do 1305 l = 1,  iback + 1
               backspace datin
1305           continue
            iback = 0
c           zero the trace array for the next panel
            DO 1320 i = 1,  ntrace
               do 1310 j = 1, ntpts
                  trace(j,i) = 0.
1310              continue
1320           continue

            go to 1000
         end if

         if(eof.or.isout.eq.irecdn) then
c           done reading events
c           close file, stop
            close(unit=datout)
            stop
         else
c           next shot
            isout = isout + 1
c           backspace once to pick up first record of new shot
            backspace datin
            iback = 0
            go to 900 
         end if
         

      end if



      if(recsrt) then

c        RECEIVER SORT
c        The number of traces is not known in advance when performing a
c        receiver sort - it depends on the recording geometry.  The
c        number of traces will be the number of different shots for which 
c        the receiver was active (and recorded some events).

c        Output one trace at a time
         ntrace = 1

         do 2000 irout = irout1,  iroutn

c        initialize the current record number to a non-existent record
         irprev = -1
c        irprev will become positive when a target event has been read
1500     read(datin,end=1750) irecd,itrace,irecv,time,
     :   xdummy,idummy,zdummy,
     :   amp,phase,caustc,head

         if(irecd.lt.irecd1) go to 1500
         if(irecd.gt.irecdn) go to 1750
         if(irecv.ne.irout)  go to 1500

         if(irecd.ne.irprev) then
c           this is a new shot (may be the first)
            if(irprev.gt.0) then
c              output the trace from the previous shot
c              write(stderr,'(1x,a,i4)')'receiver sort-record ',irprev
               write(datout)(trace(j,ntrace),j=1,ntpts)
c              zero trace for new shot
               do 1710 j = 1, ntpts
                  trace(j,ntrace) = 0.
1710              continue
            end if
            irprev = irecd
         end if

         amp = amp*fudge
         call addwav(phase,caustc,nptwav,nptby2,time,ntpts,trace,
     :   amp,wimag,wreal,whead,pi,ntrace,srate,deltat,head,resamp,
     :   maxtr,maxtpt)


         go to 1500

1750     continue
c        output last trace for this receiver sort
         if(irprev.gt.0) then
c           write(stderr,'(1x,a,i4)')'receiver sort-record ',irprev
            write(stderr,'(1x,a,i4)')'Done receiver ',irout
            write(datout)(trace(j,ntrace),j=1,ntpts)
         end if

c        done outputting sort for this receiver
c        zero trace for new receiver
         do 1810 j = 1, ntpts
            trace(j,ntrace) = 0.
1810        continue

         rewind datin
2000     continue

c        close file, stop
         close(unit=datout)
         stop
        

      end if


      STOP
      END

c-----------------------------------------------------
